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HISTORY  OF  RANDOM  NUMBER  GENERATORS 

Random  numbers  have  been  traditionally  used  in  simulation 
studies  but  have  also  proven  useful  for  solving  problems  in  numerical 
analysis.  The  Monte  Carlo  method  is  the  name  given  to  any  such  cal- 
culations involving  random  numbers .  Though  the  term  "Monte  Carlo" 

was  apparently  first  used  by  Von  Neumann  in  19^-6,  "the  idea  of  using 

2 
random  numbers  in  sampling  experiments  can  be  traced  back  to  Student 

in  1908.  Then,  and  for  some  time  later,  random  numbers  were  generated 

by  haphazard  techniques.  The  first  random  numbers  were  obtained  by 

drawing  cards  from  a  well-shuffled  deck,  throwing  coins,  or  rolling 

dice.  In  an  effort  to  reduce  the  human  bias  which  was  evident  in  such 

methods,  Tippett  collected  numbers  from  census  reports  and  published  a 

table  of  ^0,000  digits  in  1925.  Kendall  and  Babington-Smith  tested 

Tippett' s  numbers  and  found  them  inadequate  when  the  same  set  of  numbers 

was  used  over  and  over  again  in  a  sampling  experiment.  In  order  to 

supplement  Tippett' s  tables  they  published  a  new  series  of  100,000 

k 
random  numbers  produced  by  reading  digits  from  a  spinning  disk  when  it 

was  illuminated  by  a  flash  lamp.   Then  in  the  late  19^0 's  computers 

began  to  be  used  to  solve  Monte  Carlo  problems  and  so  a  need  arose  for 

a  larger  supply  of  random  numbers  of  high  quality.  The  Rand  Corporation 

filled  this  need  with  a  table  of  a  million  digits  generated  by  monitoring 

a  random  frequency  pulse  source.  This  table  was  intended  for  use  by 

punched  card  machines.  However,  computers  developed  in  sophistication 

and  it  became  impractical  for  a  high-speed  electronic  device  to  read 

random  numbers  from  storage.   It  was  more  desirable  to  have  the  computer 


2 

manufacture  its  own  numbers  as  needed  by  a  simple  arithmetic  process . 

Of  course,  numbers  produced  in  such  a  deterministic  fashion  are  not 

"truly"  random  and  so  arithmetically  generated  numbers  that  manage  to 

pass  statistical  tests  for  randomness  are  called  "pseudo"  random. 

Arithmetic  methods  are  based  on  a  recurrence 

relation  x.  ,  =  f(x.,...,x.   ),  i-n  =  0,  n  =  0,1,...  in  which  each  new 
l+l    N  i'   '  i-n"        '      '  ' 

number  is  generated  from  the  previous  ones  in  such  a  way  that  the  numbers 
appear  to  be  drawn  at  random  from  the  finite  population  of  all  the  numbers 
the  computer  can  produce.  Because  of  the  nature  of  this  recurrence 
relation,  whenever  a  set  of  n+1  numbers  is  generated  that  was  previously 
generated,  the  numbers  produced  from  then  on  will  be  identical  with  the 
numbers  produced  after  the  previous  occurrence  of  the  set  of  n+1  numbers. 
Thus  there  is  a  sequence  in  which  the  numbers  are  formed  that  will  repeat 
continuously  and  the  length  of  this  sequence  is  called  the  period. 
Starting  values  x~,  ...,x  are  needed  to  initialize  the  recurrence  relation 
and  herein  lies  one  of  the  advantages  of  arithmetic  generators.  When  an 
arithmetic  generator  is  used,  completely  new  sets  of  numbers  can  be  pro- 
duced by  changing  the  starting  values,  but  when  a  table  of  random  numbers 
is  used  the  same  digits  must  be  read  again  and  again  for  each  new  problem. 
Arithmetic  generators  are  superior  to  previously-mentioned  methods  in 
other  ways  too.  The  arithmetic  methods  are  faster,  the  generated  numbers 
are  reproducible  so  that  calculations  can  be  identically  repeated,  and 
the  numbers  can  be  analyzed  for  their  theoretical  properties  like  length 
of  period.  It  is  no  wonder  then  that  numerous  arithmetic  generators  have 
been  proposed  and  used  since  Von  Neumann  first  introduced  the  mid-square 
method  in  19^6.  In  order  to  select  the  most  appropriate  one  of  these 
methods  three  properties  of  the  generated  numbers  should  be  considered: 


3 

their  randomness,  the  length  of  the  period,  and  the  generation  time 
needed . 

In  the  mid- square  method  a  random  number  is  produced  "by- 
taking  the  middle  n  digits  of  the  square  of  the  previous  n  digit  number. 
This  new  random  number  is  then  squared  and  the  process  continues.   One 
drawback  of  this  method  is  that  a  number  may  occur  that  will  be  endlessly 
repeated.  For  example,  suppose  x  =  3500*  Then 

j?  =   12,250,000  ;  xn+1  =  2500 


Xn+1  =  6>25°>000  >     \+2   =  25°° 


Other  drawbacks  are  that  the  method  is  slow  and  the  period  is  short. 

For  these  reasons,  the  mid-square  method  was  superseded  by  the 
congruential  method.   The  congruential  method  is  based  on  the  relation 
x.  ,  =  ax.+c(mod  m)  which  means  that  ax.+c  is  to  be  divided  by  m  and  x.  n 

X+X      X  X  X+X 

n 
set  equal  to  the  remainder.  Lehmer  first  suggested  this  method  in  19^-9 

with  c=0  in  which  case  it  is  called  the  multiplicative  congruential  method. 

Q 

The  mixed  congruential  method  with  c/^0  was  later  introduced  by  Rotenberg. 
The  modulus  m,  the  starting  value  xn,  the  multiplier  a  and  the  constant  c 
are  chosen  so  as  to  provide  a  long  period  and  good  statistical  behavior. 
The  problem  of  choosing  constants  is  discussed  in  section  two.   Since 
x.  , ,  i=0,l, ...  is  equal  to  the  remainder  of  ax.+c  upon  division  by  the 
modulus  m,  each  generated  number  lies  somewhere  in  the  range  (0,1, . . . ,m-l) 
so  the  period  of  a  generator  is  less  than  or  equal  to  m.   In  this  regard 
the  mixed  method  has  an  apparent  advantage  over  the  multiplicative  method 
because  for  a  suitable  choice  of  constants  the  full  period  m  can  be 
achieved  whereas  the  full  period  can  never  be  achieved  for  the  multiplicative 


9  10 

method.  Unfortunately,  extensive  testing  '   has  shown  that  the 

multiplicative  method  generally  behaves  better  statistically  than 
the  mixed  method.  The  poor  statistical  behavior  of  the  mixed  method 
led  Hull  and  Dobell   to  conclude  that  multiplicative  methods  were  to 
be  preferred. 

In  an  effort  to  speed  up  congruent ial  methods  even  more,  the 
multiplication  operation  was  replaced  by  an  addition  so  that  the  re- 
currence relation  became  x.  ,  =   x.+x.   (mod  m) .  This  relation  is  known 

l+l    i  i-nN 

as  the  additive  method  and  when  n  =  1,  x  =  0  and  x,  =  1  one  has  the 

Fibonacci  generator.  Although  the  Fibonacci  generator  is  faster  and  has 

a  longer  period  than  the  multiplicative  method,  testing  has  shown  that 

it  gives  poorer  statistical  results.  When  n  was  increased  the  generator 

produced  numbers  with  better  statistical  properties  but  the  program 

required  indexing  so  the  speed  advantage  of  the  additive  method  was 
12 


lost, 


13 
MacLaren  and  Marsaglia   found  mixed  congruential  generators 


to  have  poor  statistical  qualities  when  used  in  Monte  Carlo  calculations 
and  so  proposed  two  improved  congruential  generators.  One  of  the  genera- 
tors uses  a  stored  table  of  random  numbers.  When  a  number  is  used  it  is 
replaced  by  a  scrambled  version  of  itself  provided  by  the  relation 
x.  =  ax.+c(mod  m) .  This  method  eliminates  the  main  problem  of  stored 
tables  of  numbers,  that  of  exhausting  the  list,  but  requires  some 
inconvenient  tape -handling.  The  second  generator  is  designed  to  improve 
the  distribution  of  n-tuples  of  random  numbers  x,,...,x  .  A  list  of 
128  numbers  u_,...,u  pn  is  produced  by  the  relation  u .  _  =   a1u_+c..(mo&  m) , 

A  second  set  of  numbers  v, ,...,v  where  v.  ,  =  a0v.+c0(mod  m)  is  also 

1'        '   n       l+l    2  l  2. 

needed.  The  first  seven  bits  of  v  are  used  to  select  the  k'th  random 


number,  x  ,  from  the  list  of  the  u's.  This  position  is  then  refilled 
with  v.  , .  Although  the  numbers  produced  by  this  generator  were  shown 

J£+_L 

to  have  "better  statistical  properties  than  numbers  produced  "by  the  mixed 
or  multiplicative  methods,  this  generator  has  other  problems.  Its 
theoretical  properties  are  unknown  and  the  generation  time  is  about  twice 
that  of  the  mixed  or  multiplicative  methods  because  two  numbers  must  he 
produced  to  generate  every  one  random  number. 

In  conclusion,  it  appears  that  there  is  no  generator  that 
provides  random  numbers  with  the  fastest  generation  time,  the  longest 
period,  and  the  best  statistical  properties.   The  generator  which 
appears  best  to  combine  these  three  properties,  though,  is  the  multi- 
plicative congruential  generator. 
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THEORETICAL  PROPERTIES  OF  THE 
MULTIPLICATIVE  CONGRUENTIAL  METHOD 


The  congruence  relation  for  the  multiplicative  method, 
x.  ..  =  ax. (mod  m),  requires  the  programmer  to  define  three  parameters: 
the  modulus  m,  the  starting  value  xn,  and  the  multiplier  a.  If  chosen 
wisely,  these  parameters  will  provide  speed,  a  long  period  and  good 
statistical  properties.  The  following  paragraphs  describe  some  rules 
by  which  the  parameters  should  be  chosen  if  these  goals  are  to  be 

attained. 

r 
The  modulus  m  is  usually  chosen  to  be  2  where  r  is  a  positive 

integer  less  than  or  equal  to  the  number  of  bits  in  a  computer  word. 
This  choice  provides  a  speed  advantage  for  suppose  that 

ax.  =  b   2r+n+...+b  2r+b  n2r"1+...+bn21+b^  ,  b  =  0  or  1. 
l    r+n         r    r-1         1    0  ' 

Now  x.  ,  =  ax.  (mod  m)  means  that  ax.  is  to  be  divided  by  m  and  x.  n  set 

l+l     ix     '  l  J  l+l 

x 
equal  to  the  remainder,  but  if  m  =  2  then  the  division  is  eliminated 

because  reduction  modulo  m  can  be  achieved  by  merely  setting  to  zero  all 

but  the  low  order  r  bits  of  ax. • 

l 

The  starting  value  xn  should  be  chosen  relatively  prime  to 

2 
the  modulus  m  or  the  period  may  be  reduced.   The  proof  of  this  statement 

requires  the  following  fact  from  number  theory.  Let  (d,m)  =  g  represent 

the  fact  that  the  greatest  common  divisor  of  d  and  m  is  g.  [Now  if 

(d,m)  =  1,  that  is  d  and  m  are  relatively  prime,  and  dx  =  dy(mod  m)  then 
x  a  y(mod  m)  if  (d,m)  =  g   1  and  dx  =  dy(mod  m)  then  x  =  y(mod  m/g)]. 
Suppose  (xn,m)  =  1  and  the  period  is  n. 


8 

x  =  ax  (mod  m)  = 

xp  =  ax.,  (mod  m)  =  a  x  (mod  m) 

• 

x  =  a  x_(mod  m) 
n     0 

Since  x  is  the  n+l'st  number  in  the  sequence  it  is  equal 

to  x  .   So  a  xQ  =  x  (mod  m)  or  a  =  l(mod  m) .  Now  suppose  (x„,m)  =  g  /  1 

and  the  period  is  n 

x  =  a  x_(mod  m) 
n      0 

a  xQ  =  xQ(mod  m) 
a  =  l(mod  m/g) 

The  period  may  be  shortened  because  the  modulus  has  been  reduced  by  the 
factor  g.  The  proof  of  this  statement  hinges  on  this  fact:  If  n  is  the 
period  modulo  m  then  it  is  the  period  modulo  m/g,  but  the  reverse  is  not 
always  true.   Suppose  a  =  l(mod  m),  then  a  =X  •  m+1  where  X   is  an  integer 
or  a  =  X.'g»m/g+l  so  a  =  l(mod  m/g).  Now  suppose  a  =  l(mod  m/g),  then 
a  =  \,»m/g+l  and  a  =  l(mod  m)  if  and  only  if  x/g  is  an  integer.  This 
fact  proves  that  the  period  when  (xn,m)  4  1  will  be  the  same  or  less  than 
the  period  when  (x„,m)  =  1.   If  the  modulus  m  is  chosen  to  be  2  ,  the 
requirement  on  x  reduces  to  choosing  x  such  that  it  is  any  odd  positive 
integer. 

The  multiplier  a  should  be  chosen  so  that  it  provides  random 
numbers  with  the  longest  period  and  the  best  statistical  properties. 
Greenberger  showed  that  the  longest  period  is  obtained  when  a  =  8k-3 

where  k  is  any  positive  integer.  The  following  results  from  number 

k 

theory  are  needed  to  explain  this  rule . 

(l)  The  Euler  Phi  function  cp(m)  is  defined  as  the  number  of  positive 
integers  less  than  m  and  relatively  prime  to  m. 


(2)  Euler's  theorem  states  that  if  (a,m)  =  1  then  q)(m)  is  the 
largest  integer  such  that  a  ^  '  a  l(mod  m) . 

(3)  If  the  least  positive  integer  n  such  that  a  =  l(mod  m)  is  such 
that  n  =  cp(m)  then  a  is  called  a  primitive  root  of  m. 

(k)     If  n  <  cp(m)  then  n  is  a  divisor  of  cp(m). 

As  was  shown  in  the  preceding  paragraph  the  period  is  the  smallest  posi- 
tive integer  n  such  that  a  =  l(mod  m)  where  we  will  assume  that  the 
modulus  m  is  2  .  We  would  like  the  smallest  n  to  he  as  large  as  possible. 

■v*  ^» 

Now  if  a  is  a  primitive  root  of  2  then  n  =  cp(2  )  by  (3)  and  if  in 
addition  a  is  odd  then  this  n  will  he  maximum  by  (2).   However,  the  only 
moduli  which  have  primitive  roots  are  h,   2»p  ,   and  p  for  p  and  odd  prime 
so  since  2  cannot  have  a  primitive  root  n  <  cp(2  ).  Now  cp(2  )  =  2    by 
(l)  and  since  a  is  not  a  primitive  root  n  must  be  of  the  form  2  for 
s  <;  r-2  by  (k) .  The  theory  requires  that  a  be  an  odd  positive  integer 
and  all  odd  positive  integers  can  be  written  in  the  form  8k-l  or  8k-3 
for  k  a  non-negative  integer.  We  would  like  to  know  which  form  of  a 
will  give  the  maximum  period. 

Suppose  that  a  is  8k-l.  The  period  is  2  for  s  ^  r-2  and  to 
determine  the  period  we  want  to  find  the  minimum  s  satisfying 

(8k^l)2  a  l(mod  2r) 
or  (8k-l)   -1  a  O(mod  2r) 

This  means  that  (8k-l)   -1  is  evenly  divisible  by  2r. 
Expanding  by  the  binomial  theorem, 

(tl+8k)2S-l  =  l2S.8k  +  2S(2S-jH8k)2  ±  ... 


=  (2s.8k) 


-1+(2S-1)  8k  +  . . . 
2 
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All  the  terms  in  square  "brackets  that  include  an  8k  will 'be  even  as 
long  as  one  of  the  denominators  does  not  cancel  out  a  numerator.  The 
fact  that  this  can  never  happen  is  seen  by  setting  n  =  2  and  rewriting 
the  j+lst  term.  We  have 

(  (n-l)(n-2)t"-  (n-J+l)  1    .  (8k)J 

The  part  in  curly  brackets  is  the  binomial  coefficient  I  n-lj  and  every 
binomial  coefficient  is  an  integer.  Therefore  every  term  in  square 


brackets  except  the  first  term  will  be  even  so  the  contents  of  the  square 

r 
brackets  is  odd.   Since  2  cannot  evenly  divide  an  odd  number  it  must  be 

true  that  2S*8k  =  0(mod  2r)   or  2S+3'k  =  0(mod  2r).  The  smallest  s 

satisfying  this  relation  depends  on  the  choice  of  k  and  is  such  that 

s  ^  r-3«  So  the  maximum  period  for  a  multiplier  of  the  form  8k-l  is 


2r"3. 


Now  suppose  instead  that  a  is  8k-3«  We  want  to  find  the  smallest 


2s 
s  such  that  (8k-3)   =  l(mod  2  ).  Again  expand  by  the  binomial  theorem 

(-3+8k)2  =  32S  +  2S'8k.32S-1  +  2s(2s-l)(8k)2.32S-2  +  ... 

2! 


2s    s 
=3   +  2s  8k 


132S-1+  (2s-l)(8k)  32S"2 


Using  the  same  reasoning  as  for  the  8k-l  case,  all  terms  in  square 
brackets  except  the  first  will  be  even  so  the  contents  of  the  square 
brackets  is  odd. 


Then  (-3+8k)2  =  2S+3.k-[odd]  +  3' 


2S        2  2s 
But  3   =  (-1+2  )   and  expanding  this  by  the  binomial  theorem 
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32S   =  1-2S.22   +  2S(2S-D(22)2 

21 


1-2S+2 


1-(2S-1)(22)  . 


The  contents   of  the  square  "brackets   is  again  odd. 

=  l-2S+2[odd] 

+  2S  r 

Now  (8k-3)   -Is  O(mod  2  ) 

g 

But  (8k-3)2   -  l=-2S+2[odd]  +  2S+3-k.[odd] 
=2S+2(-[odd]  +  2k- [odd]) 

=2S+2[odd] 

r  s+2 

Since  2  cannot  evenly  divide  an  odd  term,  it  must  divide  the  2   "  term. 

The  smallest  s  satisfying  this  is  s  =  r-2.  Therefore  the  period  for  a 

+      r-2 
multiplier  of  the  form  8k-3  is  2   .   This  proves  that  the  longest  period 

for  the  multiplicative  congruential  method  is  obtained  -when  the  multiplier 

a  is  of  the  form  8k-3. 

Unfortunately,  one  of  the  properties  of  the  multiplicative  method 

is  that  only  the  most  significant  bit  has  the  maximum  period  and  the 

length  of  period  decreases  with  the  significance  of  the  bit.  This  is  ex- 

r  4- 

plained  in  the  following  way.  Let  x    =  ax  (mod  2  ).   If  a  =  8k-3  then 

r-2 
the  period  of  the  x's  is  of  length  2   .   Suppose 

ax  =  ...b  n2r+1+b  2r+b  n2r"1+...+b  21  +  h2° 
n       r+1      r    r-1         1      0 

then  ax  (mod  2r)  =  b  12r'1+...+bn21+h2° 
nv      '  r-1         1    0 
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orx  ,  =  "b  ..2    +  y  where  y  =  ax  (mod  2   ). 
n+1    r-1      Jn      Jn     nv        ' 

r-3  r-2 

Now  the  sequence  of  y's  has  the  period  2   ,  hut  y  =  h  ?2    +  z 

r-2\  T-k- 

where  z  =  ax  (mod  2   )  so  the  sequence  of  z's  has  the  period  2 

r-i-1 
Clearly  then,  "b   .  has  period  2     .  Finally,  since  the  multiplier 

and  the  starting  value  are  odd,  each  generated  number  will  be  odd  and 

b  will  always  be  1. 

Although  a  multiplier  of  the  form  8k-3  will  provide  random 

numbers  with  the  longest  period,  these  numbers  may  not  have  good 

statistical  properties.  Unfortunately,  there  are  no  simple  rules  to 

ensure  the  latter  and  more  work  needs  to  be  done  in  this  area.  From  a 

theoretical  standpoint  the  multiplier  should  not  be  close  to  a  simple 

5 

rational  multiple  of  the  modulus  m  or  serial  correlation  will  result. 

If  a  =  m/k  for  k  an  integer  then  x.  ,  =  (m/k)x.(mod  m)  so  x.  ..  =  x./k. 

i 

Contrary  to  the  recommendation  of  some  authors  the  multiplier  should  not 
be  close  to  the  square  root  of  the  modulus  m  either,  because  this  choice 
results  in  strong  correlation  of  lag  2.  If  a  =  vm/k  for  k  an  integer  then 

x.  ,  =   (vm/k)x.  (mod  m) 

x.  2  =  (vm/k)  x.(mod  m)  =   (m/k  )x.(mod  m) 

so  x.  „  =  x./k 
i+2    i' 

Of  course  no  final  selection  of  the  multiplier  can  be  approved  until  a 

careful  check  has  been  made  on  the  generated  numbers.  For  this  purpose 

a  series  of  appropriate  statistical  tests  is  discussed  in  the  next  section. 
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STATISTICAL  TESTS 

The  methods  discussed  so  far  are  designed  to  produce  uniformly- 
distributed  random  integers.  In  addition,  a  set  of  uniformly  distributed 
random  floating  point  numbers  on  the  unit  interval  (0,1)  can  be  produced 
by  interpreting  each  binary  integer  as  a  binary  fraction  and  making  the 
appropriate  conversion.  Once  a  supply  of  uniformly  distributed  random 
numbers  on  the  unit  interval  is  available  random  numbers  "with  other 
statistical  distributions  can  he  easily  obtained.   The  hypothesis  to  be 
tested,  then,  is  that  the  generated  numbers  are  independent  and  uniformly 
distributed.  For  each  statistical  test  the  expected  results  if  this 
hypothesis  is  true  are  compared  with  the  actual  results  by  a  Chi-Square 
test. 

In  order  to  use  the  Chi-Square  statistic  the  generated  numbers 
x,,...,x  must  be  separated  into  distinct  classes  where  the  kinds  of 
classes  vary  with  each  statistical  test.   Suppose  there  are  k  cells 
r,,...,r  in  which  the  generated  numbers  can  fall  and  let  x.er.  represent 
the  fact  that  x  .  is  contained  in  cell  r. .  Let  p.  be  the  expected  proba- 
bility that  x.er.  where  j  =  1, ...,N$  i  =  1, ...k  and  p.  is  independent 
of  j.  Let  f .  be  the  actual  frequency  of  generated  numbers  in  cell  i. 

Then  X^  =  Z  (f . -N«p.)/N«p.  has  a  Chi-Square  distribution  with  k-1  degrees 

of  freedom  if  the  hypothesis  of  uniform  randomness  is  true. 

These  results  can  be  extended  to  more  than  one  dimension.  Let 

x.,...,x.   ..  be  an  m  sequence  of  generated  numbers.  As  before,  let  there 
y  j+m-l 

be  k  cells  rn,...,r,  in  which  an  individual  number  can  fall.  There  are 
1      k 
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then  k  cells  s..,...,s  ,  in  -which  an  m  sequence  can  fall  where 


(x.,...,x.   -,)es.  =  (r  ,r  ,  ...,r  )  represents  the  fact  that  the  first 
v  y        '   j+m-1'  1   v  x  y'     z 

member  of  the  sequence  fell  in  cell  x,  the  second  member  fell  in  cell  y 
and  the  id  'th  member  fell  in  cell  z.  Let  p.  be  the  expected  probability 
that  x.,...,x.   -.es.,(j  =  1, ...,N-m+l;  i  =  1,  ...,k  )  when  the  hypothesis 
is  true.  Let  f .  be  the  actual  number  of  m  sequences  of  generated  numbers 
which  are  in  class  i.  When  the  m  sequences  are  taken  consecutively, 

i.e.,  (x1,...,xm),  (xm+1> '  -  >x2x) >    (x2m+l'  *  *  *  ^  >  "  '   then 
m 

2  2 

X   =Z   (f .  -[N/m]«p.)  /[N/m]  «p.  has  a  Chi-Square  distribution.   However, 

.111       -i   m    |        J-  _1_  _L 

when  the  m  sequences  are  taken  successively,  i.e.,  (x, ,...,x  ), 

x2' ' ' '  >   x  -i )  >    (xq^  •  •  •  )t  •  •  •  "then  the  probability  that  x  . , . . .  ,x  .   _  €s  . 

km 
2  2/ 

is  not  independent  of  j  so  X   =  E  (f . -[N-m+l].p.)7[N-m+l]#p.  does  not 

i=l  1  X 

2 
have  a  Chi-Square  distribution.  For  this  case  ,  Good  shows  that 

X  -"*   .  has  asymptotically  a  Chi-Square  distribution  with  km-km~ 

2 
degrees  of  freedom.  By  defining  X   =  0,  the  case  for  m  =  1  holds 

as  before. 

2 
Once  the  Chi-Square  statistic  X  has  been  calculated  the  next 

r>  oo 

step  is  to  evaluate  /  p  f(x)dx  where 

JX 

n/2-1  -x/2 

f(x)  =  ¥-rr+n — is  'the  Chi-Square  probability  distribution  for  n 

2n/dT(n/2) 


degrees  of  freedom.  This  integral  denotes  the  probability  that  a 

2 
Chi-Square  random  variable  will  be  greater  than  X  ,  or  equivalently, 

2 
denotes  the  area  under  the  Chi-Square  curve  to  the  right  of  X  .  For 

degrees  of  freedom  n  ^  30  the  probability  can  be  found  in  a  table  of  the 
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Chi -Square  distribution.  For  degrees  of  freedom  n  >  30  the  quantity 

v2X  is  approximately  normally  distributed  about  mean  ^2n-l  with  unit 

3 
variance.   A  significance  level  of  .05  or  .01  is  usually  selected  and 

if  the  probability  is  less  than  this  significance  level  then  the  hy- 
pothesis of  uniform  randomness  is  rejected. 

Various  statistical  procedures  that  employ  such  a  test 
have  been  suggested  for  testing  the  randomness  of  the  generated  numbers. 
In  1938  Kendall  and  Babington-Smith  proposed  four  widely  used  tests : 
the  frequency  test,  serial  correlation  test,  poker  test,  and  gap  test. 
Since  then  a  vast  number  of  other  tests  have  also  been  proposed.  Many 
authors  feel  that  exhaustive  testing  of  a  generator  is  impractical,  that 
it  is  better  to  select  a  generator  that  satisfies  the  more  straight- 
forward tests  like  the  frequency  and  serial  tests.  However,  it  is  the 
opinion  of  this  author  that  the  greater  the  number  of  tests  that  are 
passed,  the  more  confidence  can  be  placed  in  the  generated  numbers. 

One  of  the  most  important  statistical  tests,  the  frequency 
test  shows  how  close  the  generated  numbers  are  to  being  uniformly 
distributed.   In  this  test  the  unit  interval  is  divided  into  100  equal 
cells  and  the  frequency  of  numbers  in  each  cell  f.,  i  =  1,...,100  is 
calculated.   If  the  numbers  are  uniformly  distributed  then  each  digit 
should  occur  an  equal  number  of  times  so  the  probability  that  a  generated 
number  is  in  any  cell  is  l/lOO.  For  a  sample  of  N  numbers  the  expected 
frequency  of  generated  numbers  in  each  cell  is  N/100.   The  hypothesis 
of  uniform  distribution  can  be  tested  by  comparing  the  calculated  fre- 
quencies with  the  expected  frequencies  by  a  Chi-Square  test.  The  sta- 

2  10°     /    2 

tistic  X  ,  =•  (IOO/N)  Z  (f.-N/lOO)  has  a  Chi-Square  distribution  with 

1  i=l  1 
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99  degrees  of  freedom.   Note  that  this  procedure  tests  the  uniform 
distribution  of  the  first  two  digits  of  every  generated  number.  It  is 
a  much  stricter  test  than  the  usual  procedure  of  dividing  the  unit 
interval  into  10  subintervals  and  thereby  only  testing  the  frequency 
of  the  first  digit  of  every  number. 

A  frequency  test  can  also  be  performed  on  the  octal  digits  of 
each  random  integer.   Suppose  there  are  m  possible  octal  digit  positions 
in  a  random  integer  and  let  f . .  be  the  frequency  of  the  i'th  octal 
digit  in  the  j'th  digit  position,  i  =  0, ...,7j  J  =  1, ...,m.   If  the 
integers  are  uniformly  distributed  then  each  octal  digit  should  occur  an 
equal  number  of  times  so  the  expected  frequency  of  each  octal  digit  in 
each  digit  position  is  N/8.  A  Chi-Square  test  can  be  used  to  compare  the 
expected  frequencies  with  the  calculated  frequencies.  The  statistic 


tions  j  =  1, ...,m  has  a  Chi-Square  distribution  with  7  degrees  of  free- 


X  .  =(8/n)E   (f . .-N/8)  calculated  for  each  of  the  m  octal  digit  posi- 

.J   \    '   -;_ n    1.1 


dom. 

The  serial  correlation  test  determines  whether  any  digit  tends 

to  be  followed  by  any  other  digit.  This  test  is  defined  by  dividing 

the  unit  interval  into  10  equal  subintervals  and  letting  f . .  be  the  fre- 

ij 

quency  of  generated  numbers  in  the  i'th  interval  which  are  followed  by  a 
generated  number  in  the  j'th  interval.  Suppose  a  sample  of  N  generated 
numbers  is  tested.  If  the  numbers  are  truly  random  then  no  digit  will 
tend  to  be  followed  by  any  other  digit  and  the  frequency  of  random  num- 
bers in  each  cell  will  be  the  same  N/lOO.   The  expected  frequencies  can 
be  compared  with  the  calculated  frequencies  by  a  Chi-Square  test.  Let 

2  10  9  o         I         \10  ,        2 

X  P  =.100/N  Z   (f,  ,-N/lOO)  and  tn  =  lO/N'Z   (f  -N/10)   vhere  X2, 
^   '     i,j=l  1J  1       i,j=l  1 
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corresponds  to  a  frequency  test  on  the  random  numbers  -with  10  divisions 

2 
of  the  unit  interval.  Then  X  -X   has  asymptotically  a  Chi-Square 

distribution  with  90  degrees  of  freedom.  A  visual  test  for  serial 
correlation  is  derived  by  noting  that  the  above  procedure  is  equivalent 
to  dividing  the  unit  square  into  100  equal  cells.  If  successive  pairs 
of  numbers  are  plotted  as  points  in  the  unit  square  then  there  should  be 
an  even  scattering  of  points  in  all  areas  of  the  square.  Any  preponderance 
of  points  in  an  area  represents  serial  correlation. 

The  mutual  independence  of  the  octal  digits  in  each  random 
integer  can  be  checked  by  a  poker  test.  To  perform  this  test,  the  first 
five  octal  digits  of  each  integer  are  treated  as  a  poker  hand  without  re- 
gard to  suit.   Each  hand  is  tallied  as  either  a  bust  (abcde),  one  pair 
(aabcd),  two  pair  (aabbc),  three  of  a  kind  (aaabc),  full  house  (aaabb), 
four  of  a  kind  ( aaaab ) }   or  five  of  a  kind  ( aaaaa ) .  For  instance ,  if  the 
first  five  octal  digits  are  62521  this  hand  is  tallied  as  one  pair.  The 
expected  frequencies  of  poker  hands  assuming  uniform  distribution  and 
mutual  independence  can  be  computed  with  the  use  of  the  multinomial  dis- 
tribution. The  multinomial  distribution  is  defined  as  follows.  If  there 
are  n  independent  trials  with  k  outcomes  per  trial  and  0.,i=l, ...,k  is 

the  probability  of  the  i'th  outcome  and  a.,i=l, ...,k  is  the  frequency  of 
the  i'th  outcome  in  n  trials  then  the  joint  distribution  of  a  , ...,a  is: 

J.       K. 

f(an,...,a,  )  =  n!       ft. 1  . . .  fl  k   <"W   :  *i 

1     k   T~* T^   1      k 

1       k 

For  the  poker  test  n  =  5,  k  =  8,  6.    =   l/8  i=l, . . . ,8 

p(5  of  a  kind)  =  8      y.  (1/8)5  [(l/8)°  •••  (l/8)°] 

(aaaaa)         51(0! • • «0!) 
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p(4  of  a  kind) 
(aaaab) 

p(Full  House)  : 
(aaabb) 


p(3  of  a  kind) 
(aaabc) 

p(2  pair) 
( aabbc ) 

p(l  pair) 
(aabcd) 

p(bust) 
(abede) 


8-7' 
8.7- 


8  A 
8  ft 


5!  (l/S^d/S)1 
TUT 

51   (l/8)3(l/8)< 
312! 


5\       (l/8)3(l/8)1(l/8)1 
3!1!1I 

_5I_  (l/S^Cl/S^d/S)1 
2I2I1. 


^    (l/8)2(l/8)1(l/8)1(l/8)1 
2111111! 


51  (1/8)5 

1! 


For  a  sample  of  random  integers  of  size  N  the  expected  frequencies  of  the 
7  poker  hands  are  N»p.  ,i=l, . . .  ,7*  The  actual  frequencies  f.  can  be 
compared  with  the  expected  frequencies  e.  by  a  Chi-Square  test.  The 


7 


■f.-e.    N2 


statistic  X     =  Z     (  i     i   )     has  a  Chi-Square  distribution  with  6  degrees 

i=l  e~. 

l 

of  freedom.   In  the  interpretation  of  this  test  a  positive  correlation 
in  the  octal  digits  'would  result  in  too  many  good  hands  while  a  negative 
correlation  would  result  in  too  many  bad  hands. 

A  run  test  shows  how  the  generated  numbers  oscillate  among 
themselves.  To  describe  this  test  suppose  we  replace  the  sample  of  N 


random  numbers  x  , . . .  ,x  by  an  N-l  sequence  S=s  , . . .  ,  s    where 

s.=0  if  x.  <x.  .  and  s.  =  1  if  x.  >x.,n.  Runs  can  be  ins ide  runs  or  end 
i      i  ^  l+l      i        i    l+l 

runs  depending  on  whether  they  occur  inside  or  at  the  two  ends  of  the 
N-sequence.  A  subsequence  of  k  zeros  bounded  by  ones  at  each  end  forms 
an  inside  run  up  of  length  k.  A  subsequence  of  k  ones  bounded  by  zeros 
at  each  end  forms  an  inside  run  down  of  length  k.  End  runs  have  to  be 
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"bounded  only  on  one  side.  For  example,  the  sequence  001110  begins  with 
a  run  up  of  length  2  and  is  followed  by  a  run  down  of  length  3*  The  run 
test  involves  counting  the  number  of  runs  of  different  lengths  and  com- 
paring them  by  a  Chi-Square  test  with  expected  results.   In  the  interpre- 
tation of  this  test  a  high  frequency  of  long  runs  indicates  non-randomness 
of  the  generated  numbers .  The  expected  number  of  runs  of  length  k  is : 

2N.k  +3k+l  -2,  k3+3k2-k-*i- 
(k+3)!     (k+3): 

This  is  derived  in  the  following  manner.   Let  r  be  the  number  of  runs 

of  length  k. 

Let  r,  .  =/  1  if  there  is  a  run  of  length  k  starting  in  the  i'th  position 

of  the  S  sequence 

0  otherwise 

Then  the  expected  number  of  runs  of  length  k  is: 

N  N 

E(r  )  =  E(Z  rki).  Now  E(r^)  =  0,  i  >  N-k  so  E(Z  r   )  = 
i=l  i=l 

N-k-1 

E(rn  ,  )  +  Z  E(r,  .  )  +  E(r.  ,T  ,  ).  For  2  <  i  <  N-k-1  then 
x  kl   ,  p  v  ki'      k,N-k'        v   x 

E(r,  .)  =  PClO1^)  +  P(01k0)  =  2-P(lOkl) 

=  2-P(x._1  >  x.  <C  x.+1  <C  •  •  •  £  x1+krl  £  x.+k  >  x.+k+1) 

1 

Since  0  <  x .  <  1,  then  P(x  . )  =   /  dx . 
J  J     Jri   J 


P(X.  ,  <  X.)  a  /   JdX.  _ 


X 

0 
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r1 

P(x.   ..    >  x.)   =      /     dx.   , 
J"1         J  Jx       J"1 

j 


So  E(r.  , . )  =  2   /  "  '  f'  I  'Xi+k.  •  •  /  Xl+1  f   **±^3x±~  '^i+k-l^i+k^i+k+l 

i+k+1  x 

=  2-k2+3k+l 
(k+3)! 

Nov  E(r      )  =  P(Okl)   +  P(lkO)  =  2-P(0kl) 

=  2*2{x±  ^  x2  ^   . . .  ^  xk  <  xk+1  >  xk+2) 

=  2    /      /  /  /        . . .    /        dx_,  . . .  dx,    ,  dx,  dx,    ,  dx. 


1  k-1     k     k+1     k+2 

0     x,  .  n     „  0  0 


k+2     0 


=  2.  k+1 

Xk+2T: 


But  since  E(r,  ,  )   =  E(rv  w  t^   "then 

N 
E(rk)   =  E(Z_  rk.)   =  S.E^)   +   (N-k-2)E(rkj_) 


=  2N  k2+3k+l  -  2.k3+3k2-k-+ 
(k+3)!  (k+3)! 


Another  type  of  run  test  counts  the  number  of  runs  above  and 
below  the  mean.  To  describe  this  test  suppose  we  replace  the  sample  of  N 

random  numbers  x, , . . .  tx«,  by  an  N  sequence  S  =  s, , . .  .  ,s.T  where  s .  =  0 
1'    'TI  IN        l 

if  x.  ^  l/2  and  s.  =  1  if  x.  >  l/2.  The  number  of  runs  up  and  down  of 
different  lengths  are  counted  as  described  in  the  previous  paragraphs 

and  compared  with  expected  results  by  a  Chi-Square  test.  The  expected 

-k-1 
number  of  runs  of  length  k  is  (N-k+3)2    .  This  is  calculated  in  the 
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following  manner.  Let  r  "be  the  number  of  runs  of  length  k. 

Let  r..  =(l  if  there  is  a  run  of  length  k  starting  in  the  i'th  position 

of  the  S-sequence. 

0  otherwise 

Then  the  expected  number  of  runs  of  length  k  is : 

N  N 

E(r,  )  =  E(Z  r     ).     Now  E(r   . )  =  0,   i  >  N-k+1  so  E(Z  r,  . )  = 
K  1=1  1CL  Kx  i=l  K± 

N-k 

E(rkl}   +  &>*(**)  +  E(rk,N-k+l)'     *or  2  <C  i  <;  N-k  then 


E 


(r.  .)   =  P(lOkl)  +  P(01k0)  =  2.p(lOkl) 


ki 


=  2.P(x._1  >  1/2,   x..   <C  1/2,...,^  <  1/2,   x.+k  >  1/2) 


But  P(x.  ^  1/2)   =  1/2,   P(x.  >  1/2)   =  1/2 
J  J 

So  E(rk±)   =  2-P(2~k"2)   =  2'k_1 

Now  ECr^)   =  F(o\)   +  P(lk0)   =  2-P(0kl) 

=  2»P(x1  ^  1/2, ...,xk  <C  1/2,   xk+1  >  1/2)   =  2.(2-k_1) 
=  2"k 


But   since  ECr^)   =  E(rk  N_k+1)  "then 

N 
E(rk)   =  E(Z  rfcL)   =  S.ECr^)   +   (N-k-l)E(r]d) 

=   (N-k+3)2"k_1 


A  test  which  is  similar  to  the  run  test  is  the  gap  test.  If  a 
maximum  in  a  sequence  of  generated  numbers  is  a  number  that  is  surrounded 
by  two  smaller  numbers  and  a  minimum  is  a  number  surrounded  "by  two  larger 
numbers  then  a  gap  is  defined  as  the  distance  between  two  successive 
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maxima  or  minima.  In  determining  the  length  of  a  gap  the  two  maxima  or 

minima  are  included  as  -well  as  the  numbers  "between  them  so  the  smallest 

gap  is  of  length  3.   The  gap  test  is  performed  by  counting  the  number  of 

gaps  of  different  lengths  and  comparing  them  with  expected  results  by  a 

Chi-Square  test.  For  a  sample  of  generated  numbers  of  size  N  the  expected 

number  of  gaps  of  length  r  is : 

r-1 
3«N-2 -  r-2  .  This  formula  is  derived  very  clearly  by  Kermack- 

rl    r+2 

McKendrick  and  it  will  not  be  discussed  here. 

7 
MacLaren  and  Marsaglia  were  dissatisfied  with  the  afore- 
mentioned statistical  tests  because  they  thought  the  tests  were  irrelevant 
to  actual  uses  of  random  numbers.   One  of  the  most  common  ways  of  using  a 
sequence  of  random  numbers  is  to  sample  consecutive  n-tuples  and  so 
MacLaren  and  Marsaglia  devised  frequency  tests  on  the  range  of  n-tuples 
and  on  the  range  of  maximums  and  minimums  of  n-tuples.   If 

¥.  =  max(x.  _,..., x.   )  i  =  0, . ...N-n  where  the  x. .  are  independent  and 
1      v  l+l'    '  i+n7      '   '  l's 

uniformly  distributed  on  the  unit  interval  (0,1)  then  W.  is  not  uniformly 
distributed  because : 

P(W.  <  a)  =  P(x.  .  <  a)'--  P(x.   <  a) 

v  1    '  v  l+l    '  v  1+n    ' 

pa  pa 

=   /   dx.  .  . . .  /   dx. 
J0    i+l    JQ         i+n 

n 

=  a  , 

whereas  uniformity  would  require  P(W.  <  a)  =  a.   However, 

P(W.  <  a)  =  P(W.   <  a  )  =  a  so  ¥.   is  theoretically  uniformly  dis- 
11  l 

tributed  on  (0,1).   In  a  sample  of  N  generated  numbers  there  are 

N/n  n-tuples  and  the  observed  W.  ,  i  =  1, ...,N/n  can  be  tested  for  a 

uniform  distribution  by  the  frequency  test  as  described  earlier.   In  a 


2k 


similar  manner  the  minimum  of  n-tuples  can  be  tested.  In  this  case,  if 

W.  =  min(x.  .,,..., x.   )  then 

1      x  i+l    '  i+n' 

P(W.  <  a)  =  1-P(W.  >  a)  =1-  [p(x.+1  >a)     P(xi+n  >  aij 


=  !■ 


dx 


a 


i+l 


J 


■1 


dx. 


a 


l+n 


,n 


n  . 
is 


=  l-(l-a) 

Now  P(Wj,  <  a)  =  P  l-fl-W^11  <  l-(l-a)n   =  l-(l-a)n  so  l-(l-W±) 
theoretically  uniformly  distributed.  A  frequency  test  can  again  be 
employed  to  determine  how  close  the  actual  results  are  to  being  uniformly 
distributed. 

Unfortunately,  a  random  number  generator  which  has  passed  all 
of  these  standard  statistical  tests  may  still  fail  when  a  particular 
application  is  sensitive  to  a  property  of  the  numbers  that  was  not 
tested.   In  the  end,  the  most  relevant  test  of  a  random  number  generator 
for  a  particular  application  is  to  use  the  generator  on  a  similar  problem 
for  which  the  answer  is  known. 
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EXPERIMENTAL  RESULTS 

The  multiplicative  method  of  random  number  generation  has  been 
generally  accepted  as  a  good  one  because  its  theoretical  properties  are 

known,  it  is  fast  and  easy  to  program,  it  has  a  long  period,  and>  most 

1  2 
important,  it  produces  numbers  with  good  statistical  properties.  ' 

Unfortunately,  published  results  of  the  application  of  the  multiplicative 
method  so  far  have  been  for  computers  with  word  lengths  of  36  bits  or 
larger.   There  is  some  speculation  whether  the  multiplicative  method  will 
be  an  adequate  method  of  random  number  generation  for  the  IBM  360.  A  prob- 
lem inherent  within  the  method  concerns  the  selection  of  a  multiplier 
which  will  produce  numbers  with  good  statistical  properties.  From  a 
theoretical  standpoint  it  is  known  that  the  multiplier 

1)  should  be  of  the  form  8k-3  for  k  a  non-negative  integer  in 
order  to  ensure  the  maximum  length  of  period 

2)  should  not  be  close  to  a  simple  rational  multiple  of  the 
modulus  or  serial  correlation  will  result 

3)  should  not  be  close  to  a  simple  rational  multiple  of  the  square 
root  of  the  modulus  or  correlation  of  lag  2  will  result 

but,  aside  from  these  rules,  little  is  known  about  why  some  multipliers 
produce  numbers  with  better  statistical  properties  than  others.  This  prob- 
lem is  magnified  by  the  3^0 's  small  word  length  of  32  bits  and  it  is 
believed  that  random  number  generation  on  the  3^0  by  the  multiplicative 
method  will  be  extremely  sensitive  to  choice  of  a  multiplier. 

In  an  attempt  to  learn  more  about  random  number  generation  on 
the  IBM  360  the  multiplicative  method  was  programmed  for  the  360.  A  total 
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of  1,500  different  multipliers  was  generated  randomly  and  used  to 
generate  sequences  of  pseudorandom  numbers.   Each  sequence  was  statis- 
tically tested  and  the  test  results  were  compared.   Integer  arithmetic 
was  used  to  generate  fixed  point  numbers  and  floating  point  numbers 
on  the  unit  interval  were  obtained  by  setting  to  zero  the  first  eight  bits 
of  each  integer  number  and  inserting  the  appropriate  exponent.   In  order 
that  the  bit  patterns  for  the  integer  random  numbers  and  floating  point 
random  numbers  be  identical  for  the  statistical  tests  the  integers  were 

stored  as  twenty  four  bit  results  with  the  first  eight  bits  of  each 

2k 

word  set  to  zero.  The  modulus  of  the  method  thus  was  reduced  to  2 

and  since  the  starting  value  used  was  odd  and  all  multipliers  were  of  the 

form  8k-3  for  k  a  non-negative  integer  the  period  of  all  the  generators 

22 
tested  was  the  maximum  possible  of  2   .  To  test  the  1,500  random 

multipliers  the  decimal  number  2173  was  arbitrarily  selected  as  the 
starting  value  and  for  each  multiplier  a  sequence  of  5>000  or  10,000 
numbers  was  generated.   The  frequency  test  and  serial  correlation  test 
were  performed  on  each  sequence  of  generated  numbers  and  for  both  tests 
a  Chi-Square  statistic  and  the  associated  Chi-Square  probability  (i.e., 
the  probability  that  if  the  numbers  are  random  a  Chi-Square  variable  will 
exceed  this  test  statistic)  were  calculated. 

Of  the  1,500  multipliers  tested  the  ten  which  had  the  largest 
Chi-Square  probabilities  were  selected  for  further  testing.   These  ten 
multipliers  were  compared  in  three  different  ways  with  the  nine  multi- 
pliers which  gave  the  lowest  Chi-Square  probabilities  and  for  which  the 

hypothesis  of  uniform  randomness  was  rejected.  First,  since  a  multiplier 

2k 

could  be  as  large  as  2  -1  it  could  have  as  many  as  eight  decimal  digits 

and  it  was  thought  that  larger  decimal  numbers  were  probably  better 


28 

multipliers  than  smaller  ones.  However,  some  of  the  five -digit  numbers 

were  "better  multipliers  than  the  eight -digit  numbers  so  this  idea  had 

3 
to  be  rejected.   Coveyou  stated  that  the  multiplier  should  not  have  a 

small  number  of  l's  and  it  was  conjectured  by  this  author  that  the  "best 
multipliers  probably  had  an  equal  number  of  O's  and  l's  in  the  binary 
representation.  After  tallying  O's  and  l's  it  appeared  that  this  idea 
also  had  to  be  rejected.  It  was  then  suggested  that  "by  looking  at  the 
octal  digits  of  the  multipliers  it  might  be  observed  that  certain  octal 
digits  were  present  or  were  lacking  in  particular  digit  positions  of  the 
ten  acceptable  multipliers.  This  hypothesis  again  proved  false  and  so 
the  search  for  finding  rules  to  predict  when  a  multiplier  would  be  a 
good  one  was  concluded. 

Next  it  was  desired  to  learn  more  about  the  kinds  of  generators 
that  were  produced  by  using  the  ten  selected  multipliers  as  defined  in 
the  previous  paragraph.   Passing  tests  like  the  frequency  and  serial 
correlation  tests  is  necessary  but  certainly  not  sufficient  for  the 
acceptance  of  a  generator  so  each  multiplier  was  subjected  to  an  exten- 
sive set  of  tests.  In  all,  sixteen  tests  were  performed  on  the  first 
10,000  numbers  generated  with  starting  value  2173*  The  tests  are 
described  as  follows : 

TEST  1  -  Frequency  test  on  the  generated  numbers  with  100  divisions  of 

the  unit  interval. 
TEST  2  -  Serial  correlation  test.  Successive  2-tuples  of  numbers  were 

taken  as  points  in  the  unit  square.  The  unit  square  was  divided 

into  100  equal  cells  and  the  frequency  of  9>999  pairs  was 

tallied. 
TEST  3  -  Successive  3 -tuples  of  numbers  were  taken  as  points  in  the  unit 
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TEST  7 


eube.  The  unit  cube  was  divided  into  1,000  equal  cells  and 
the  frequency  of  9>998  triples  was  tallied. 

TEST  k   -  Consecutive,  non-overlapping,  2-tuples  of  numbers  were  taken 

as  points  in  the  unit  square.  The  unit  square  was  divided  into 
100  equal  cells  and  the  frequency  of  5>000  pairs  was  tallied. 
\V  TEST  5  -  Consecutive,  non-overlapping,  3-tuples  of  numbers  were  taken  as 

points  in  the  unit  cube.  The  unit  cube  was  divided  into  1,000 
equal  cells  and  the  frequency  of  3>333  triples  was  tallied. 

TEST  6  -  Consecutive,  non- overlapping,  4-tuples  of  numbers  were  taken 
as  points  in  the  unit  cube  in  four  space.  The  unit  cube  in 
four  space  was  divided  into  10,000  equal  cells  and  the  fre- 
quency of  2,500  4-tuples  was  tallied. 

Frequency  test  on  the  eight  octal  digit  positions  of  each 
generated  integer  number. 
Poker  test. 

Test  of  runs  up  and  down. 
Test  of  runs  above  and  below  the  mean. 

The  maximum  number  in  consecutive,  non- overlapping,  2-tuples  of 
numbers  was  calculated.  The  unit  interval  was  divided  into  100 
equal  cells  and  the  frequency  of  the  5>000  maximums  was  tallied. 

TEST  12  -  The  maximum  number  in  consecutive,  non -overlapping,  3-tuples 

of  numbers  was  calculated.  The  unit  interval  was  divided  into 
100  equal  cells  and  the  frequency  of  the  3>333  maximums  was 
tallied. 

TEST  13  -  The  maximum  number  in  consecutive,  non- overlapping,  ^-tuples 

of  numbers  was  calculated.  The  unit  interval  was  divided  into 
100  equal  cells;  frequency  of  the  2,500  minimums  was  tallied. 


<] 


,vV> 
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TEST  8 


TEST  10 
TEST  11 
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For  each  test  the  calculated  results  for  the  generated  sequence 
of  numbers  were  compared  with  results  that  would  have  "been  expected  from 
a  truly  random  sequence  of  numbers  "by  a  Chi-Square  test.   The  results  of 
the  Chi-Square  test  are  given  in  TABLE  1.   Each  column  of  the  table 
represents  a  particular  multiplier,  shown  in  decimal  at  the  head  of  the 
column,  and  each  row  represents  a  particular  test.  The  numbers  listed 
represent  the  Chi-Square  probability  or  percentage  probability  of  exceed- 
ing the  calculated  Chi-Square  statistic.   The  appearance  of  e  in  the  table 
means  that  the  Chi-Square  probability  was  less  than  .05  and  that  the 
hypothesis  of  uniform  randomness  was  rejected.  The  eight  probabilities 
listed  for  TEST  7  are  the  results  of  the  frequency  test  performed 
separately  on  each  of  the  eight  octal  digit  positions  of  the  2k   bit  random 
integers.  The  first  number  is  for  bits  22-2^  (counting  beginning  from  the 
right  end  of  the  word),  the  second  number  bits  19-21,  and  the  eighth 
number  bits  1-3 •   It  is  obvious  that  the  Chi-Square  probability  for  the 
eighth  digit  position,  bits  1-3,  will  be  extremely  low  since  the  least 
significant  bit  must  always  be  a  1  because  the  generated  numbers  are  odd. 
The  two  tests  on  the  digits,  tests  7  and  8,  are  not  as  important  as  the 
other  tests  because  most  applications  require  numbers  rather  than  digits. 
In  this  regard,  the  manner  in  which  TABLE  1  is  viewed  depends  largely  on 
the  particular  application  for  which  a  generator  is  needed.   If  an 
application  requires  randomness  of  consecutive  2 -tuples  then  the  results 
of  tests  k,   11,  and  1^  indicate  that  multiplier  9  should  not  be  chosen. 
If,  instead,  good  randomness  of  consecutive  3-tuples  is  needed  then  tests 
5,  12,  and  15  indicate  that  any  of  the  multipliers  is  adequate.   Tests  6, 
13,  and  l6  show  that  for  randomness  of  consecutive  ^-tuples  multiplier  10 
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TEST  1 
TEST  2 
TEST  3 
TEST  k 
TEST  5 
TEST  6 
TEST  7 


TEST  8 
TEST  9 
TEST  10 
TEST  11 
TEST  12 
TEST  13 
TEST  Ik 
TEST  15 
TEST  16 


TABLE  1 

-1- 

13651723 

-2- 
12316023 

-3- 

50515 

-k- 
29355 

-5- 
5^27 

.9* 

.<* 

.98 

•  97 

•  99 

.90 

.92 

.98 

•97 

.89 

.93 

•  57 

•  98 

Al 

.76 

-99 

.68 

.87 

•  99 

.57 

•  79 

•  55 

.kQ 

•  53 

.k2 

.80 

.27 

M 

•  33 

.83 

•  95 
.1£ 

•  51 

•  99 

•  99 

•  99 

•  99 

.85 
.76 
.16 

•97 
•  99 
•99 
•99 

.76 

•  55 
.k2 

•  99 
•99 

•  99 
-99 

.81 

•  77 

•  98 
•99 

•  99 
-99 
-99 

.69 
.50 

•  58 
-99 
-99 

•  99 

•  99 

6 

G 

e 

e 

€ 

.08 

.68 

.07 

-19 

.27 

e 

.07 

•35 

.56 

.11 

.26 

.ko 

•93 

.76 

.90 

•97 

•96 

•97 

.83 

AO 

.9^ 

•  57 

.62 

•  75 

-95 

.65 

.86 

.96 

.29 

.18 

•  53 

•  05 

.98 

•  37 

.92 

•27 

•  30 

-9^> 

.76 

•  99 

.90 

.13 

.88 

.72 

.7* 
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TABLE  1  (Continued) 


TEST  1 
TEST  2 
TEST  3 
TEST  k 
TEST  5 
TEST  6 
TEST  7 


TEST  8 
TEST  9 
TEST  10 
TEST  11 
TEST  12 
TEST  13 
TEST  Ik 
TEST  15 
TEST  16 


-6- 
223227 

-7- 
220915 

-8- 

220615 

-9- 
3739387 

-10- 
1^59251 

•  95 

•  98 

.<* 

.96 

•96 

•  93 

.89 

•  93 

.88 

.98 

•  52 

•  39 

.88 

.81 

.60 

.kk 

.89 

•  51 

.72 

•91 

.kd 

.66 

•  93 

.92 

.3^ 

.76 

.90 

•  97 

.62 

•95 

•  59 

•  55 
.65 

•  99 

•  99 

•  99 

•  99 

.83 

•  55 
.89 

•  93 

•  99 

•  99 

•  99 

.72 

•  3^ 

.66 

•  99 

•  99 

•  99 
-99 

.84 

.63 
.09 

•  99 

•  99 

•  99 

•  99 

.5^ 
•93 

•  38 

•  99 
-99 
-99 
-99 

e 

€ 

£ 

€ 

e 

•15 

.kk 

Al 

G 

.92 

.85 

•32 

.8lt- 

e 

.20 

.68 

.Ik 

•  27 

.28 

•89 

•  52 

.71 

•  70 

e 

.ko 

.60 

.kz 

-9^ 

.22 

•  57 

.54 

.k6 

.66 

.06 

.kl 

•  19 

•  98 

.06 

.85 

•15 

•  70 

.84 

.73 

.kd 

.21 

.1+8 

.70 

.35 

.17 

e 
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should  not  "be  chosen.  Some  applications  require  uniform  numbers  that 
oscillate  randomly  among  themselves.  This  property  of  the  generated  num- 
bers is  tested  for  by  tests  9  and  10  and  the  results  of  these  tests  show 
that  multipliers  1  and  9  should  not  be  used  for  this  application. 

Although  the  ten  carefully  chosen  multipliers  of  TABLE  1 
gave  good  results  it  must  be  remembered  that  the  10,000  numbers  tested 
in  each  case  comprise  a  very  small  sample  of  the  ^,19^,30^-  possible 
numbers  in  the  period  of  the  generator.   Consequently,  further  tests  were 
performed  using  different  sample  sizes  and  different  starting  values. 
TABLE  2  shows  the  effect  of  a  change  in  starting  value.  Each  row  repre- 
sents a  particular  decimal  starting  value  shown  at  the  left  of  each  row 
and  each  column  represents  a  particular  decimal  multiplier  shown  at  the 
head  of  the  column.  The  numbers  listed  in  the  table  are  the  Chi-Square 
probabilities  for  the  frequency  test  performed  on  a  sample  of  10,000 
numbers  generated  with  each  particular  starting  value .  TABLE  3  shows  the 
effect  of  a  change  in  sample  size.   Each  row  represents  a  particular  size 
sample  shown  at  the  left  of  the  row  and  each  column  represents  a  particu- 
lar multiplier  shown  in  decimal  at  the  head  of  the  column.   The  numbers 
listed  in  the  table  are  the  Chi-Square  probabilities  for  the  frequency 
test  performed  on  particular  size  samples  generated  with  a  starting  value 
of  2173.  As  the  statistics  in  TABLES  2  and  3  show,  when  one  repeats  the 
same  test  on  different  samples  one  should  get  a  range  of  Chi-Square 
values  with  the  corresponding  Chi-Square  probabilities  distributed  uni- 
formly. Since  the  interpretation  of  the  statistics  does  not  vary 
drastically  with  a  change  in  sample  size  or  starting  value  there  is  no 
reason  to  believe  that  the  original  tests  on  samples  of  10,000  numbers 
with  starting  value  2173  were  biased. 
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This  "battery  of  tests  was  next  used  to  learn  more  about  the 
IBM  random  number  generator  for  the  360,  RANDU.  There  had  "been  some 
hesitation  to  use  RANDU  because  its  statistical  properties  were  not  known. 
The  same  sixteen  tests  were  applied  to  the  first  10,000  numbers  generated 
by  RANDU  with  starting  value  2173  and  the  Chi-Square  probabilities  are 
listed  in  TABLE  4  under  multiplier  11.  The  results  were  theoretically 

predictable.  The  generator  RANDU  uses  the  multiplicative  method  with  the 

31 
modulus  2   and  the  multiplier  65539*  Note,  however,  that  65539  is 

l6 

2   +3*  IBM  made  a  mistake  in  choosing  the  multiplier  so  close  to  the 

square  root  of  the  modulus  and  this  accounts  for  the  extremely  high 
correlation  among  3-tuples  as  shown  by  the  results  of  TEST  3» 

One  of  the  most  recent  developments  in  the  random  number  genera- 

4 
tion  field  is  a  paper  by  Marsaglia  in  which  he  reports  that  all  multi- 
plicative generators  have  a  defect  in  that  when  successive  n-tuples  of 
numbers  (x_ , . . . ,x  ),  (xp, . . . ,x   ), . . .  are  plotted  in  the  unit  n-cube 
all  the  points  will  fall  into  a  small  number  of  parallel  hyperplanes.  This 
defect  is  not  too  serious  since  most  users  will  not  select  random  numbers 
in  this  manner  anyway.  If  an  application  does  require  randomness  of 
successive  n-tuples,  however,  some  improved  methods  of  generation  have 
been  suggested.  One  of  these  methods,  suggested  earlier  by  MacLaren 
and  Marsaglia,  involves  the  mixing  of  two  separate  multiplicative  genera- 
tors. This  method  was  programmed  for  the  360  and  produced  random  numbers 
x.,...,x  in  the  following  manner.  A  list  of  128  numbers  u_,  ...,11.-0  vas 

24 
produced  by  the  relation  u.  1  =  a  u.  (mod  2  ).   Another  relation 

24 
v.   =  a  v. (mod  2  )  was  used  to  generate  a  second  set  of  numbers 

v,  ,...,v  .  Bits  18-24  of  vl were  used  to  select  the  k'th  random  number, 
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TEST  1 
TEST  2 
TEST  3 
TEST  k 
TEST  5 
TEST  6 
TEST  7 


TEST  8 
TEST  9 
TEST  10 
TEST  11 
TEST  12 
TEST  13 
TEST  Ik 
TEST  15 
TEST  16 


TABLE  4 

/2#i>J>4 

MWi 

-11- 

65539 

-12- 

220651 
13651723 

.6? 

.9k 

.38 

.58 

\ 

.86 

.85 

•  90 

•  31 

.k6 

.07 

•  91 

•  32 

•  93 
.20 

-99 
■  99 

•  99 

•  99 

.96 
.61 
•  51 
.99 
-99 
-99 
-99 

•  Ik 

.09 

.21 

.20 

.26 

.50 

.64 

•  70 

•  71 

•97 

.81 

.21 

•  97 

•3^ 

•  71 

.6k 

•  37 

.68 
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x  ,  from  the  list  of  the  u's  and  this  position  was  then  refilled  with 


v   .   The  two  multipliers  selected  for  a  and  a  were  multipliers  8  and  1 
respectively.  The  sixteen  tests  were  performed  on  the  first  10,000 
numbers  generated  with  starting  value  u  =  2173  and  v  =  2173  and  the 
Chi-Square  probabilities  are  listed  in  TABLE  k-   under  the  heading  multi- 
plier 12.  Though  this  method  eliminates  the  inherent  non- randomness  of 
successive  n-tuples  in  the  multiplicative  method  it  requires  128  storage 
locations  and  is  slower  than  the  multiplicative  method.  The  time  for 
generation  of  one  random  number  for  multipliers  1-10  and  RANDU  was  .0022 
hundredths  of  a  second,  but  the  generation  time  for  the  Marsaglia  mixing 
process  was  .0031  hundredths  of  a  second.  With  the  speed  and  storage 
capacities  of  modern-day  computers,  however,  these  two  problems  should  not 
be  serious  for  most  users . 

In  conclusion,  a  battery  of  programs  was  written  to  provide 
stringent  statistical  tests  for  a  random  number  generator.  A  group 
of  random  number  generators  was  devised,  each  one  meeting  the  needs  of  a 
different  type  of  application.   It  is  hoped  that  the  user  will  consult 
TABLES  1-k   and  then  select  the  generator  which  best  satisfies  his 
requirements .   In  the  likely  event  of  a  further  breakthrough  in  the  field 
of  random  number  generation  the  battery  of  statistical  programs  will 
undoubtedly  be  of  assistance  in  testing  future  random  number  generators. 
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FOOTNOTES 

1.  Coveyou  and  MacPherson,  J.  Assoc.  Comp.  Mach.  Ik. 

2.  Hull  and  Dobell,  J.  Assoc.  Comp.  Mach.  11. 
3«   Coveyou  and  MacPherson,  lk. 

k.     Marsaglia,  Proc.  N.  A.  S.  6l.     d ****•»*  ■ 

5.  MacLaren  and  Marsaglia,  J.  Assoc.  Comp.  Mach.  12. 
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APPENDIX 


^3 


SUBROUTINE  NAME  ***  RANB 

IDENTIFICATION  ***  Random  number  generator 

PURPOSE  ***  Generates  a  sequence  of  uniformly  distributed  random 

floating  point  numbers  on  the  unit  interval  (0,1)  and  uniformly 

22 
distributed  random  integers  on  the  interval  (0,2   ).   The  period 

22 
of  the  generated  numbers  is  2 

OTHER  SUBROUTINES  USED  ***  None 

USAGE  ***  REALM  RN(NX) 

INTEGER  *4  IRN(NX) 


CALL  RANB  (NSTART,N,RN,IRN) 
Where  the  parameters  in  the  calling  sequence  are : 

NSTART  -  The  starting  value.   It  should  be  an  odd  positive  integer. 
N  -  The  number  of  random  numbers  to  be  generated,  NX  >  N. 
RN  -  The  one  dimensional  array  in  which  the  N  floating  point  random 
numbers  are  stored  by  RANB.   Its  contents  may  be  arbitrary 
upon  entry. 
IRN  -  The  one  dimensional  array  in  which  the  N  integer  random  numbers 

are  stored  by  RANB.   Its  contents  may  be  arbitrary  upon  entry. 
Note :  When  a  particular  multiplier  is  desired,  the  card  labelled 

MULT  in  the  program  should  be  changed  to  the  correct  multiplier. 
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SAVhS     1HI-    STARTING    VALUt- 


7HKMS    lllll     THt-    t-jKSl     H    HITS    HI-    THH     I'MTHGhK 
SI  i  IK  I- S     TMh     iNfhKPR     Ii\i    Ihi-    AKKAY     IKN 
CIinVPR'I      I  li    PLUA'f  I>\'fi    H C 1 1  ivi  I 

f;«-  I  S     IHf-    KAiVUtiH    iMUmMfcR     I'M    A     hLTINU    PI     RK; 
NliKf-  AL  1  Zl-S     THt-    KAMDilW    MHKiHhK 

GETS    THfc    NWF'tritn*    rtA(.K     I'M    GbNL.    KPG.    b 
PUTS     KAlMDOivi     IMUmHpK     Jim     I  HP     AKKAY     K  IM 
VII  IS     THp     I'M'IbOK     I  l\i    KH;.     4 


KhSlUWhS     THt-     SIAKlI'v(;     \/A|.U(- 


ChAkACTHK IS1 1L 
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SUBROUTINE  ***  UNINUM 

IDENTIFICATION  ***  Frequency  test  for  random  numbers 

PURPOSE  "**■*  Determines  how  close  a  set  of  random  numbers  are  to  being 
uniformly  distributed  by  dividing  the  unit  interval  into  100  subinter- 
vals,  counting  the  frequency  of  random  numbers  in  each  subinterval 
and  comparing  with  expected  results  by  a  Chi-Square  test. 

OTHER  SUBROUTINES  USED  ***  PROLIM 

USAGE  ***  REAlM  RN  (NX) 

CALL  UNINUM  (N,RN,PR0B) 

Where  the  parameters  in  the  calling  sequence  are: 

N  -  The  number  of  random  numbers  to  be  tested,  NX  ^.  N. 

RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 
tested. 

PROB  -  The  outputed  value  of  the  Chi-Square  probability.   It  is  a 

number  between  0  and  1  which  indicates  the  probability  that  a 
Chi-Square  variable  will  exceed  the  calculated  test  statistic. 


k6 


S  l  •  •-»  W  I  H  II   IWl-     NHlPVlIN,  (  |\.t  KM,  PWIIH  ) 

IhSI     '»»-    O'vlHIKM'fY    hF    NUmH^RS    WITH    100    DlVlSIIMS    UF    1 HK    UNIT     INTERVAL 

H(=4L*4    HM(w|fHIV|lfil) 
?iM|  R;^k*<.    (>LLI  100) 

niv(i  )«o.() 

dii  ^  T=2.mi 

oi w(  I  )=mv(  T-i  )  +  .oi 

on    h    l=l,) oo 
Cf-LL(J)=n 

III  I      H       [  =  1   ,  I>1 

DII    7    .1=1,100 

1 1-  (w>'(  I  )  .<;i-.i>]\/(.|)  .ANI),KN(  I  )  .LT.DIVJ  J  +  l  )  )f;n    TH    « 

L.'jNl  IMIh 

CHLL  (  J  )  =CHLL  ( .1  )  +  1 

|-XPMlS  =  M/]00. 
XrtJllHsM.O 

mi  w   i  =  i , loo 

XW(ji»isXNt)M+(CHLL  (  I  )-t-XM0IS)**2 

iiim(. S"=xi\i(ii"i/eXh,ii!  s 

CaLL    ijkiiL  !""<  ww,  I  hvf.Su,  PWIIR) 

MHTUKim 

(-Ml) 


SUBROUTINE  NAME  ***  FRESOC 

IDENTIFICATION  ***  Frequency  test  on  successive  or  consecutive  n-tuples 

of  random  numbers . 
PURPOSE  ***  Determines  how  close  a  set  of  n-tuples  of  random  numbers 

are  to  being  uniformly  distributed.  The  n-tuples  can  be  taken  either 

successively  or  consecutively  and  n  must  be  1,  2,  3>  or  4.  The  unit 

2 
(interval,  square,  cube  or  cube  in  4-space)  is  divided  into  (10,  10  , 

3       h- 
10  ,  or  10  )  subintervals  respectively  depending  on  the  value  of  n. 

The  frequency  of  random  numbers  in  each  sub interval  is  tallied  and 

compared  "with  expected  results  by  a  Chi-Square  test. 

USAGE  ***  REAlM  RN(NX) 

GALL  FRESOC  (%RN?NTUPLE?=IS0C=CHISQ^PR0B) 

Where  the  parameters  in  the  calling  sequence  are : 

N  -  The  number  of  random  numbers  to  be  used  in  the  test,  NX  ^N. 

RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 
tested. 

NTUPLE  -  The  size  of  the  tuple,  NTUPLE  =  1,2,3,4. 

ISOC  -  If  ISOC  =  1  then  the  N  random  numbers  are  to  be  divided  into 
N/NTUPLE-eonsecutive,  non- overlapping,  tuples  of  numbers. 
CHISQ  -  On  exit  from  FRESOC  CHISQ  will  contain  the  calculated  Chi- 
Square  statistic  for  the  test.   If  ISOC  =  1  the  input  value  for 
CHISQ  must  be  the  CHISQ  value  obtained  from  a  previous  call  to 
FRESOC  with  the  NTUPLE  value  one  less  than  the  NTUPLE  value  for 
this  call.   If  the  NTUPLE  value  for  this  call  is  1  then  the  input 
value  for  CHISQ  is  0.   If  ISOC  =  2  the  input  value  for  CHISQ  may 
be  arbitrary. 


PROB  -  The  outputed  value  of  the  Chi  Square  probability.   It  is  a 
number  between  0  and  1  which  indicates  the  probability  that  a 
Chi -Square  variable  will  exceed  the  calculated  test  statistic 
CHISQ. 
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SOBkOOl  I  .Mi-    |-wt-SOC(N,RN,l\!TOPLF,  I  SoC,OK'.Su,  PwilB  ) 

FRFOUHNCY     I>S1     UN    SUCCESSIVE    OR    CONSFCUTIVF    N    TUPLES 

KFftL*4    HN(W)  ,  UI  V(  11  ) 

Iivjl  FGF«*4    CFLL  (  10,  10,10,10)  ,LPM{4) 

l)TV(  1  )=o.o 

l)D    5    1  =  /,  1  1 

I )  T  w  (  T  )  =  i )  1  v  I  I  - 1  )  +  0  . 1 

J  =  l  _ 

K=l 

L  =  l 

mi  y   i=i,io 
no  r   j  =  i  ,  l  o 

1)11    7    K=l  ,  10 
Dll    ft    L=l,10 
CFLLI  I  ,.),  K,l.  )=0 
IFJNTUPLE.FO.  1  )GU    TO    9 
IF(NTUPLF.FO.2)G0    fo    fi 
IF(n"IUPLE.F0.3)GU    TO    7 
COM  IIiMtIF 
COMT  IiMllfc 
CON  "I  IlMllh 
CUM  T  IIMOH 

00  10     1=1,4 
LPM( I )=l 

IF(  ISOC.fcu.2  )(,n"  TO  To" 
INC«=1 

1  k  M  A  X  =  m  -  im  IUPLb+I 
Gil    TO    ?>i 
MOi»iTOP  =  lM/ivTOPLt:" 

IKMAX  =  imiimiip#nT0PLF-imTUPLF+_1 

INCR=NTUPLE 

DO  50  I«=l , IkMAX, INCR 

00  40    MT=1 .  fivTltPLF 
IRP= IH+MJ-I 

DO     30    J=l  ,  10 

IRKM  IRP)  ,GE.l>TV(  J)  ,AMD.RN(  IRP)  .LT.DI  V(  J+l  )  )G(.I    TO    40 

CONTINUE 
„LPM(NT)=J 

CFLULPMI  1  )  ,LPM(  2)  ,LPM(  3)  ,  LPM(4)  )=CtLL(LPM(l)  ,LPM(2)  ,LPM(  3)  ,LPM(4) 
1  )  +  l 

1  F  (  IS0C.bi-).2)G0    TU~~bi~ 

FXPDIS=lRMAX/(  10.**IMTUPLE) 

GU  TO  54 

FXPDIS=M0mT0P/ ( 10.*#NTUPLE) 
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s^ 


S7 


NUI'I 
=  1 
=  1 

=  1 

n  *6 

II     b 

n  5 
N 1 1 M 
F  (  M 
F  (  N 
F  (  M 
[}i\iT 


=  0.0 


0  1=1,10 
4  -1  =  1,10 
H  K  =  l  ,  10 
7    L=l,10 

=XNUM+(CELL(  It  J»K,L)-EX  PETS') 
"I  UPLE.EQ.  1  )0(i    Til    60 
riJPLE.E0.2)r,f)    TU    b9 
rUPLE.E0.3)Gf]    TM    68 

1  M  I  I  h 


**2 


5H  COMTIMUE 

^V  CUNT  INI  J  E 

60  CHN'I  I  IMliE 

IE(  IS!)C.EQ.2)GD     In    yo 

TEwiCSO=XMUM/EXPI)I  S 

MmCS"=  I  thCSO-iiiMC  S'i 

i-i  i ) I-  =  1 0  *  * N T U P L  E- 1 0**  (  N f  UPL  E- 1 

GO    T(l    HO 

HMCSn=XNl.iM/EXPI)TS 

Mf)F=]  0**n  riJPLfc 

CALL     ^41  iL  J  M  (  NDH  ,  UNCStl  ,  PROFS  ) 

R  F  T UR M 

FlMI'l 
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SUBROUTINE  NAME  ***  SERPLT 
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IDENTIFICATION  ***   Plotter  routine  for  random  number  generator 

PURPOSE  ***  Successive  2 -tuples  of  random  numbers  are  plotted  as  points 
in  the  unit  square .  Any  preponderance  of  points  in  an  area  of  the 
plot  indicates  serial  correlation.  Note  that  the  numbers  printed  along 
each  axis  should  be  divided  by  10  to  obtain  their  actual  values . 

OTHER  SUBROUTINES  USED  ***  CCP1PL,CCP5AX 

USAGE  ***  REAlM  RN(NX) 

GALL  SERPLT (N,RN) 
Where  the  parameters  in  the  calling  sequence  are : 
N  -  The  number  of  random  numbers  to  be  used  in  the  test,  NX  ^.  N. 
RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 
tested. 


C". 

40 

50 
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S  lift  ROUTINE    SERPLT  (  N,RN) 

REAL*A    RN(N) ?T(2) .EDGE 

[mTEGER*4  UP, DOWN, LAST ( 16,16) 

LGRID=16 _. _ 

C      THINK  BASK  LGRID 

C      LGRIO  MUST.BE  _POWER  OF  TWO_ 

IMJ  10  1=1, LGRID 

DO  10  J  =  l, LGRIO       

10     LAST(T,J)=0 

Mnls=N-i 
C      PUT  POINTS  IN  BOXES 

on  so  1=1, NM1 

LAB=LGR  Il)*RN~(  I  )  +  1 

LORD=LGRID*RN( 1+1 )+l 

LPRFV=LAST< LAB,LOWU) 

fh(LPKEV.E0.0)G0  TO  40 

LDIFF=I-LPREV 

R  N  (  I  )  =  R  N  (  I  )  +  L  D  I  F  F 

STORED  BACK  POINTER  IN  RN(  I  j 

LAS  f  ( LAB, LORD) =1 

CI  ImT  INUE 

HAVF  POINTS  IN  BOXES,  NOW  DECODE  AND  PLOT. 

I  ( ]  )=0.0 

T(2)=1.2b 

CALL  CCPlPL(0.b,0.t>,-3) 

CALL  CCPbAX(0.0,0.0f 21HSEC0ND  MEMBER  OF  P A  I  R  ,  2  1  ,  8 . 0 , 90 . 0 , T ) 

C/U.L  CCPbAX ( 0.0,0.0, 20HFIRST  MEMBER  OF  PA  I  R  ,  -20  ,  B  .  0,  0  .  0  ,  "I  ) 

0P  =  3 

i  ii  IWN=2 

FDGF=0.01 

DO  100  LAB=1,LGRID 

OH  99  L0RD=1,LGR ID 

I=LAST (LAB, LORD) 
60      1  F  (  1  .Hi.O  )G0  TM  99 

JX=RM( I ) 

JY  =  RN'(  1  +  1  ) 

RN( 1 )=RN{ I )-JX 

BY=RN( 1+1 )-JY 

XC=B.O*RN( I ) 

YC=B.O*FY 

XL=XC-EDGE 

Xm=XC+EDGF 

YL  =  YC-EDGF 

Y*«=YC  +  EDGE 

CALL  CCP1PL ( XL, YL, DP ) 

(.ALL  CCP1  PL  (  XM,  YW»  DOWN) 

CALL  CCP1PL ( XM, YL,OP ) 

CALL  CCP1PL ( XL , Ymf DOWN) 

In! PP=JX 

IF(  I  DIFF.  h(i.O  )GI  I  H>  99 

I  =  I-|J)IFF 

GO  id  60 
99  l".i  i.  •  I  T  •Of- 
1  00    CI  I  I  1 1  < I  IE 

RETURN 

f  (v  n 
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SUBROUTINE  NAME  ***  POKDIG 

IDENTIFICATION  ***  Poker  test  and  frequency  test  on  the  octal  digits  of 
a  set  of  random  integers. 

PURPOSE  ***  The  2k   bit  random  integers  are  divided  into  8  octal  digits  each. 
In  the  poker  test  the  first  5  octal  digits  of  each  integer  are  treated 
as  a  poker  hand  without  regard  to  suit.  The  frequencies  of  certain 
poker  hands  are  tallied  and  compared  with  expected  results  by  a 
Chi-Square  test.   In  the  frequency  test  the  frequency  of  each  octal 
digit  in  each  of  the  8  digit  positions  is  tallied  and  compared  with 
expected  results  by  a  Chi-Square  test. 

OTHER  SUBROUTINES  USED  ***  FREPOK, PROLIM 

USAGE  ***  REAlM  PR0Bl(8) 

INTEGER  *K   IRN  (NX) 

• 

CALL  POKDIG(N, IRN, IFOP, PR0B1,PR0B2) 
Where  the  parameters  in  the  calling  sequence  are : 
N  -  The  number  of  random  integers  to  be  used  in  the  test,  NX  >  N. 
IRN  -  The  one  dimensional  array  containing  the  N  random  integers  to 

be  tested. 
IFOP  -  If  IFOP-1  then  the  frequency  test  is  to  be  performed.  If 

IFOP  =  2  then  the  poker  test  is  to  be  performed. 
PR0B1  -  If  IFOP  =  1  then  PROBl(l),  I  =  1, ...,8  contains  the  Chi-Square 
probability  of  the  frequency  test  for  the  i'th  digit  position 
where  the  1st  digit  position  is  bits  22-24  (counting  from  the 
right  end  of  the  word),  the  2nd  digit  position  is  bits 
19-21,...,  the  8'th  digit  position  is  bits  1-3.  The  Chi-Square 
probability  is  a  number  between  0  and  1  which  indicates  the 


5^ 
probability  that  a  Chi -Square  variable  "will  exceed  the  calculated 
test  statistic.   If  IFOP  =  2.,   then  PROKL(l),  I  =  1,...,8  will  be 
arbitrary  upon  exit  from  POKDIG. 
PR0B2  -  If  IFOP  =  1  then  PR0B2  will  be  arbitrary  upon  exit  from 

POKDIG.   If  IPOP  =  2  then  PR0B2  contains  the  Chi-Square  probability 
for  the  poker  test.   It  is  a  number  between  0  and  1  which  indicates 
the  probability  that  a  Chi-Square  variable  will  exceed  the  calcu- 
lated test  statistic. 
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4 


14 


lb 


20 


26 


1000 


SUBROU1  IMP     POKDIG  (N,  IRN,  IF0P,PRr>Bl,PK(JB2  ) 
FREQUENCY    TEST    AMI)    POKER    TEST    UN    THE    01  GITS 
REAL*4    PXPPOK ( 7  )  , PRUB1 ( 8 ) 
INTEGERS    MCT ( H, 8) ,POK(7) , IRN(N) 


i) 


1  =  1,  H 


5 

Dfl     3     J  =  l,8 

()CT(  I  ,J)=0 

Dfl    4    1  =  1,7  

PMK  (  T  )=() 

CALL     PR  F  Pf.  IK  (  N  ,  I  RN.flLQ  J_tRJJKjLl  FOP-l: 

I F ( IFUP.E0.2 )GU    TO    20 

FXPFRE=)M/8, 

00    lb    J=1,R 

XNUM=0 , 0  . 

DU     14     1=1, H 

XNUM=XNUM+(UC I ( I , J )-PXPPRE)**2 

CHI  FRE  =  Xi\iUM/EXPERE 

CALL    PR OL I M ( 7 , C H I  PR P , P  R  U  B 1 ( J )  ) 

GO    IN     1000 

,?0i>07R*N 

i>12695*N 

lb3809#l\! 

102b39*N 

01 7090#N 

008b45*N 

000244*IM 


EXPPUK ( 1 )=0 
PXPPOK ( 2 ) =0 
PXPPUK ( 3 )=0, 

PXPPMK ( 4) =0 

PXPPUK  (  b)=0. 

PXPPMK (6)=0 

PXPPUK ( 7 )=0, 

CHI  PMK  =  0.() 

Ull    2  b    1  =  1,7 

CHIP(IK  =  CHIPl)K+(PHK(  I  )-PXPP()K(  I  )  ) 

CALL    PRULIto{6,CHIP0K,PR0B2) 

RETURN 

FNO 


2/PXPPl)K(  I  ) 
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SUBROUTINE  NAME  ***  FREPOK 

IDENTIFICATION  ***  A  subroutine  called  by  POKDIG 

PURPOSE  **"*  The  2k   bit  random  integers  are  divided  into  8  octal  digits 
each.  In  the  poker  test  the  first  5  octal  digits  of  each  integer  are 
treated  as  a  poker  hand  without  regard  to  suit.  The  frequencies  of  cer- 
tain poker  hands  are  tallied  by  FREPOK.   In  the  frequency  test  the  fre- 
quency of  each  octal  digit  in  each  of  the  8  digit  positions  is  tallied 
by  FREPOK. 

OTHER  SUBROUTINES  USED  ***  None 

USAGE  ***  INTEGER*^-  IRN(NX),0CT(8,8),P0K(7) 

• 

CALL  FREP0K(N,IRN,0CT,P0K,IF0P) 
Where  the  parameters  in  the  calling  sequence  are : 
N  -  The  number  of  random  integers  to  be  used  in  the  test,  NX  ^  N. 
IRN  -  The  one  dimensional  array  containing  the  N  random  integers  to  be 

tested. 
OCT  -  If  IFOP  =  1  this  two  dimensional  array  is  set  to  zero  on  entry 

to  FREPOK.   On  exit  from  FREPOK  OCT(l,j);I,J  =  1, ...,8  contains 

the  frequency  of  the  octal  digit  1-1  in  the  J'th  digit  position 
where  bits  22.-2h   are  the  1st  digit  position,...,  bits  1-3  are  the 

8th  digit  position.   If  IFOP  =  2  the  contents  of  this  array  are 

arbitrary  on  entry  and  exit  from  FREPOK. 
POK  -  If  IFOP  =  2  then  POK(l),  I  =  1, . . . ,7  is  set  to  zero  on  entry 

to  FREPOK.   On  exit  from  FREPOK  POK  contains  the  frequencies 

of  7  poker  hands  where : 

POK(l)  =  frequency  of  bust 

P0K(2)  =  frequency  of  one  pair 
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P0K(3)  =  frequency  of  two  pair 
POK^)  =  frequency  of  3  of  a  kind 
P0K(5)  =  frequency  of  a  full  house 
P0K(6)  =  frequency  of  k   of  a  kind 
P0K(7)  =  frequency  of  5  of  a  kind 

If  IFOP  =  1  then  the  contents  of  this  array  are  arbitrary  on  entry 
and  exit  from  FKEPOK. 
IFOP  -  If  IFOP  =  1  then  the  frequency  of  digits  is  to  be  tallied.   If  IFOP 
2  then  the  frequency  of  poker  hands  is  to  be  tallied. 
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PUTS    N     IN    REG.     b 

PUTS  .ADORES S__UF. OCT  [  1 ,1 )  IN  R EG ... ...  1 0 

PUTS  ADDRESS  UP  IRN(l)  IN  REG.  11 

PUTS  IFOP  IN  REG.  4 


THIS  INITIALIZES  THE  DIST  ARRAY  T(l  ZE 

THE  DJJjT  ARRAY.  STORES  THE  NUMBER 

OF  TIMES  THAT  EACH  OCTAL  DIGIT  UCCURR 

IN  THE  EIRST  FIVE  OCTAL  DIGITS  (IE 

A  RANDOM  NUMRER. 

PUTS  IRN( I )  IN  REG.  3 

GETS  RID  OF  INITIAL  ZEROS  OE  IRN(I) 

REG.  9  TELLS  WHICH  OCTAL  POSITION 

R  EG .  7  TELL  S  WH I CH  OCTAL  DIGIT 

CLEARS  REG.  2 

MOVES  AN  OCTAL  DIGIT  INTO  REG.  2 

HOUND  WHICH  OCTAL"  DIGIT  IT  IS 


RO 
ED 


EREUUENCY  TEST  ONLY 


ADDRESS  MODIFICATION  OE  OC T ( 1 , 1 ) 

RO IS  OCT ( I , J )  IN  RED.  2 


RES  TOR  PS  REG.  9 


MURE  OCT  AL  DIGITS  IN  IRM  I  ) 

POTS  ADDRESS  Oh  ROK(l)   IN  REG.  2 

RFG.  h    TELLS  WHETHER  3  OF  A  KIND 

REG.  7  TELLS  MOW  MANY  PAIRS 

REG.  H  USED  10  INCRFMENI  THh  DIST  ARRAY 
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SUBROUTINE  NAME  ***   RNSUAD 

IDENTIFICATION  ***  Test  of  runs  up  and  down  on  a  set  of  random  numbers 

PURPOSE  ***  If  x..,...x  are  the  set  of  random  numbers  then 

x  ,  ^  x  >x  _>•••>  x  ,  <  x  ,  , .   is  a  run  down  of  length  k, 
n-1    n    n+1        n+k    n+k+1  ' 


and  x  -.  >  x  ^  x  n<«««^x  n  >  x  ,  -   is  a  run  up  of  length  k. 
n-1    n  ^  n+1  **    ^  n+k    n+k+1  ° 


This  test  involves  tallying  runs  up  and  down  of  different  lengths 
and  comparing  with  expected  results  by  a  Chi-Square  test. 
OTHER  SUBROUTINES  USED  ***  PROLIM 


USAGE  ***  REAlM  RN(NX),EXPRUN(100) 
INTEGER^  RUNLTH(IOO) 


CALL  RNSUAD(N,RN,EXPRUN,RUNLTH,PR0B) 

Where  the  parameters  in  the  calling  sequence  are : 

N  -  The  number  of  random  numbers  to  be  tested,  NX  ^  N. 

RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 
tested. 

EXPRUN  -  Internal  array  used  by  RNSUAD,  its  contents  are  arbitrary 
on  entry  to  RNSUAD.  On  exit  from  RNSUAD  it  contains  the 
expected  frequencies  of  runs  of  lengths  1, . . .m  where  m  is 
the  maximum  run  length  found  in  the  N  random  numbers. 

RUNLTH  -  Internal  array  used  by  RNSUAD,  its  contents  are  arbitrary 

on  entry  to  RNSUAD.  On  exit  from  RNSUAD  it  contains  the  actual 
frequencies  of  runs  of  lengths  1, . . .m  where  m  is  the  maximum 
run  length  found  in  the  N  random  numbers. 
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PROB  -  The  outputed  value  of  the  Chi -Square  probability.  It  is  a  number 
between  0  and  1  which  indicates  the  probability  that  a  Chi -Square 
variable  will  exceed  the  calculated  test  statistic. 
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TO     10 


SUBROUTINE    RNSUAD(N, RN, EX PR UN , RUNLTH , PROB ) 

TEST    DF    RUNS    UP     AND    DOWN 

REAL#4    RN(N) , EXPRUN( 1 ) 

INTEJ2ER.*^_RJJNLTHJ  1) 

M4XLTH=6 

DO    2     1  =  1  ,  100 

RUNLTH ( I )=0 

LAS  "I  7=0 . 

LAST  1  =  0 

1  =  1 

TF(RM(  I  ) ,GT.RN(  T  +  l  )  )G0    TO    10 

1  =  1  +  1 

IF(  I  ,EO.N)G(J    Tl)    13 

IF(RN( I ) .LE.RN( 1+1) )G0    T 0    b 

K=  I -LAST  1-1 

L  A  S  T  7  =  I  -  1 

MSWICH=1 

GO    Tm    IS 

1=1  +  1 

I  F  (  I  ,EO.I\l)G( 

IF{RN{ I ) .GT 

K=  I  -LA  ST  7.-1 

L  A  S  T  1  =  I  -  1 

NSWICH=0 

GO     III     IS 

K=M-LAST1-1 

GO    TO    IS 

K=m-LAST7-1 

RUNLTH(K  )=RUi\iLTH(  K  )  +  1 

IF ( K.G1 .MAXLTH) mAXLTH=K 

IF( I .Em. m) GU    TO    ?3 

IF(NSWlCH.EO._0JGg     Id    5 

GO    "in    10 

U(i    ^s    I  =  1,MAXLTH 

IP3=I+3 

EACT=1.0 

DO    ?A    J=l, IP3 

FACT=FACT*J 

ISU*T**2 

EXPRON( I ) =  ?.*(N  = 

I  M AX  =  i"AXLIH 

Ull     Sf)     1  =  1,  I  MAX 

I 1= lMAX-I+1 

IE  (ROiMLTHl  I  I  )  .GE.SlGO    TO     SI 

RUNLTH (  I  1-1  ) =RUNLTH (  I  I  —  1  ) +RUNLTH (  I  I  ) 

EXPRUN(  I  I -1  )  =  EXPRUN(  I  I -1 ) +EXPRUN(  I  I  ) 

MAXLTH=MAXLT H-l 

Chi  sn  =  f).o 

on    (So    I  =  i,MAXLTH 

CHI  SO  =  CHI  S0+(  RUNLTH  (  I  ) -EX  PR  UN  (  I  )  )_**2/ EXPRUNI  I  ) 

NUF=i"iAXLTH-l 

CALL  PRUL IM( MDFtCHISQf PROB) 

R  E  T  0  R  N 

E  M I ) 


(  TSU+3.*I  +  1.  )-(  I  SO*  I +  3.':=  I  SO- I -4.  )  )  /FACT 
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SUBROUTINE  NAME  ***  RNSABM 

IDENTIFICATION  ***  Test  of  runs  above  and  below  the  mean  on  a  set  of 

random  numbers. 
PURPOSE  ***  If  x, , . . .x„  are  the  set  of  random  numbers  then 

x  _  <C  l/2,  x  _  ^  l/2 ,  . . .  ,x         ^  l/2  is  a  run  below  the  mean 

of  length  k  and  x  ,  >  l/2,  x  p  >  l/2,...,x    >  l/2  is  a  run 

above  the  mean  of  length  k.  This  test  involves  tallying  runs  above 

and  below  the  mean  of  different  lengths  and  comparing  with  expected 

results  by  a  Chi-Square  test. 

OTHER  SUBROUTINES  USED  ***  PROLIM 

USAGE  ***  REAlM  RN(N),  EXPRUN(IOO) 
INTEGER  *k   RUNLTH(IOO) 

• 

GALL  RNSABM(N,RN,EXPRUN,RUNLTH,PR0B) 

Where  the  parameters  in  the  calling  sequence  are: 

N  -  The  number  of  random  numbers  to  be  tested,  NX  >  N. 

RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 
tested. 

EXPRUN  -  Internal  array  used  by  RNSABM,  its  contents  are  arbitrary  on 

entry  to  RNSABM.  On  exit  from  RNSABM  it  contains  the  expected 
frequencies  of  runs  of  lengths  1, . .  .m  where  m  is  the  maximum 
run  length  found  in  the  N  random  numbers . 

RUNLTH  -  Internal  array  used  by  RNSABM,  its  contents  may  be  arbitrary 
on  entry  to  RNSABM.   On  exit  from  RNSABM  it  contains  the 
actual  frequencies  of  runs  of  lengths  1, ...,m  where  m  is 
the  maximum  run  length  found  in  the  N  random  numbers. 
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PROB  -  The  outputed  -value  of  the  Chi -Square  probability.   It  is  a 
number  between  0  and  1  which  indicates  the  probability  that 
a  Chi-Square  variable  will  exceed  the  calculated  test  statistic. 
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£>  A  P' 

SUBROUTINE  NAME  ***  maxmin 

IDENTIFICATION  ***  Frequency  test  on  maximums  or  minimums  of  n-tuples 
of  a  set  of  random  numbers.  If  W  is  the  maximum  of  the  n  numbers  in 
an  n-tuple  then  X  =  w  is  uniformly  distributed^  if  W  is  the  minimum 
of  the  n  numbers  in  an  n-tuple  then  Y  =  l-(l-W)  is  uniformly  dis- 
tributed. In  this  test  the  unit  interval  is  divided  into  100  sub- 
intervals  and  the  frequency  of  X's  and  Y's  is  tallied  and  compared 
against  expected  results  by  a  Chi-Square  test. 

OTHER  SUBROUTINES  USED  ***  UNINUM 

USAGE  ***  REAlM  RN(NX),MX0MN(MX) 

• 

GALL  MAXMIN(N,RN,]OTUPI£,MXOMN,IXON,PROB) 
Where  the  parameters  in  the  calling  sequence  are : 
N  -  The  number  of  random  numbers  to  be  tested,  NX  ;>  N. 
RN  -  The  one  dimensional  array  containing  the  N  random  numbers  to  be 

tested. 
NTUPLE  -  The  number  of  random  numbers  in  a  tuple.  Note  that  the  tuples 

are  taken  consecutively  from  the  N  random  numbers. 
MXOMN  -  An  internal  array  used  by  MAXMIN,  MX  ^.n/NTUPLE.   Its  contents 
are  arbitrary  on  entry  to  MAXMIN;  on  exit  from  MAXMIN  it  con- 
tains the  N/NTUPLE  X's  or  n/NTUPLE  Y's. 
IXON  -  If  IXON  =  1  then  the  maximum  of  n-tuples  is  to  be  tested.   If 

LXON  =s  2  then  the  minimum  of  n-tuples  is  to  be  tested. 
PROB  -  The  outputed  value  of  the  Chi-Square  probability.  A  number 
between  0  and  1  which  indicates  the  probability  that  a 
Chi-Square  variable  will  exceed  the  calculated  test  statistic. 
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SUBROUTINE     MAXHIN(  Nf RN, N TUPLE , MXOMN,  IXON, PROB) 
REAL*4    RN(N)  ,  MXOMIM(  1)  _. 

I  M  =  N  / 1\|  i  I J  P  L  E 

I  k  iv,  a  X  =  IN*NTyPLE-NTUPLE+l 

i\lTMl=NTUPLE-I 

1  =  1 

l)(")    40    J  =  l  ,  TRMAX,NTUPLE 

MXUMM(  I  )  =«IM(  J  )  ______  ._ ,'__ 

I E( IX0N.E0.2 )GU    TO    20 

DO     10    K=ltMThl 

IF(RM(  J  +  K)  .  GT.MXOM!M(  1  )  )  MXOMN  (  I  )  =RN  (  J  +  K  ) 

MXUMNt  I  )=MXOMN(  1  )**N.TUPLE. 

GO     TO    40 

D\\    30    K  =  l  ,  iMTMl 

IF(RM( J+K) .LT.MXOMNI  I  )  )MXOMN(  I  )=RN( J  +  K ) 

MXOMN( T ) =1 .-( 1 . -MXOMN( I ) ) **NTOPLE 
1  =  1  +  1 

CALL     UN  I  IMUM  (  lNf  MXI IMM  ,  PKUH  ) 
R  E  T  U  R N 
FND 
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SUBROUTINE  NAME  ***  PROLIM 

IDENTIFICATION  ***  Chi -Square  probability 

PURPOSE  ***  Given  a  Chi-Square  statistic  X     and  the  number  of  degrees  of 

freedom  n,  this  subroutine  calculates  the  probability  that  a  Chi- 

2 
Square  variable  will  exceed  X  .  For  n  ^  30  this  is  done  by  evaluating 

/DO  II 

f(x)dx  where  f(x)  =  x     e       is  the  Chi-Square  probability 
x2  2n/2r(n/2) 


distribution.   For  n  >  30  the  quantity  n/2X  -  \/2n-l  is  normally  distributed 

f  °°  -  2/2 

so  /  p  f (x)dx  where  f(x)  =  1   e  '  ' 


X 


n/2^ 


is  evaluated. 

OTHER  SUBROUTINES  USED  ***  None 

USAGE  ***  CALL  PROLIM(NDF  CRTSQ,PROB) 

Where  the  parameters  in  the  calling  sequence  are : 

NDF  -  The  number  of  degrees  of  freedom,  an  input  parameter. 

2 
CHISQ  -  The  Chi-Square  statistic  X  ,   an  input  parameter. 

PROB  -  The  Chi-Square  probability,  an  output  parameter.   It  is  a 

number  between  0  and  1  which  indicates  the  probability  that  a 

Chi-Square  variable  will  exceed  the  statistic  CRTSQ. 
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SI  muni  JUNE    PRML  Im(  NUF,UNCSQ,  PROB) 
K(LAL*4    X(96')fW(%)  

INC1 (  I  )  =  (?**<  (NUH-2)/2.  )  )*EXP(-T/2.  ) 

IMC2 ( T )=bXP( (-1**2) /2. ) 

01)  =    00. 0162767 A 

02)  =    -0,01627674 

03)  =  00.04Rfi].29R 

04)  =  -0.04HR129H 
Oh)  =  00. OR  129749 
06)  =  -0.08129749 
07  )  =  00. 1 1 369585 

08)  =  -0.11369585 

09)  =  00.1459737] 

10)  =  -0.14597371 

11)  =  00. 17R0968R 

12)  =  -0.178096R8 

13)  =  00.21003131 

14)  =  -0.21003131 
l  b  )  =  00.241 74315 

16)  =  -0.24174315 

17)  =  00.27319881 
in)  =  -0.27319RR1 

19)  =  00.30436494 

20)  =  -0.30436494 
21  )  =  00. 3352085  2 

22)  =  -0.33520852 

23)  =  00.36569686 

24)  =  -0.36569686 

25)  =  00.39579764 

26)  =  -0.39579764 

27)  =  00.4254789  8 
?H)  =  -0.42547898 

29)  =  00.45470942 

30)  =  -0.45470942 
3]  )  =  00.48345797 

32)  =  -0.48345797 

33)  =  00.51169417 

34)  =  -0. hi  1694  17 

35)  =  00.53938810 

36)  =  -0.53938810 

37)  =  00.36651041 

38)  =  -0.56651041 
3S  )  =  00.59  30  32  36 

X(40)  =  -0.59303236 
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